%%%%%%%%%%%%%%%%%%%ITC分析
% addpath(genpath('C:\Users\custs\workspace\matlab\congcong\lib\eeglab2023'));
addpath(genpath('C:\Users\custs\workspace\matlab\congcong\lib\fieldtrip-master'));
cd 'C:\Users\custs\workspace\matlab\congcong\data\aMCI\clean'
if ~isempty(findall(0, 'Type', 'Figure'))
    close all;
end
clear;
id={'103','106','204'};
Ns = length(id);

for subi=1:Ns
    load([id{subi} '_amciCleanData.mat']);

    %======== time frequency decomposition ==========
    cfg              = [];
    cfg.output       = 'fourier';
    cfg.pad          = 'nextpow2';
    cfg.channel      = 'all';
    cfg.method       = 'mtmconvol';
    cfg.taper        = 'hanning';       % one hanning taper
    cfg.foi          = 1:0.05:30;          % 可以选择更密的频率点，eg, 1:1:30
    cfg.t_ftimwin    = 0.5 .* ones(size(cfg.foi));
    cfg.toi          = -0.7:0.001:1.5;            % 可以选择更密的时间点, eg, -0.2:0.01:0.5
    fourier_A  = ft_freqanalysis(cfg, condA);    % for A condition
    fourier_B  = ft_freqanalysis(cfg, condB);    % for B condition

    %========== compute ITC ============
    % condition A
    itc_A           = [];
    itc_A.label     = fourier_A.label;
    itc_A.freq      = fourier_A.freq;
    itc_A.time      = fourier_A.time;
    itc_A.dimord    = 'chan_freq_time';
    F_A        = fourier_A.fourierspctrm;               % copy the Fourier spectrum
    itc_A.itc  = squeeze(abs(mean(F_A./abs(F_A),1)));   % compute ITPC

    % condition B
    itc_B      = itc_A ;
    F_B        = fourier_B.fourierspctrm;
    itc_B.itc  = squeeze(abs(mean(F_B./abs(F_B),1)));

    %======= baseline correction  =========
    %%% METHOD I: subtract baseline average
    %     cfg = [];
    %     cfg.parameter    = 'itc';
    % %     cfg.baseline     = [-0.2 0];
    % %     cfg.baselinetype = 'absolute';      %静息态不基线校正
    %     itc_A = ft_freqbaseline(cfg,itc_A);
    %     itc_B = ft_freqbaseline(cfg,itc_B);
    %
    %     itc_A.itc = atanh(itc_A.itc);     % Fisher-Z transform
    %     itc_B.itc = atanh(itc_B.itc);


    %     %%% METHOD II: shuffling, this is time consuming
    %     num_iter    = 500;            % number of iteration，可按需要改成1000
    %     Ltimes      = size(F_A,4);    % number of time points
    %     Ntrials     = size(F_A,1);    % number of trials
    %     itc_A_surrogate  = zeros(num_iter,size(F_A,2),size(F_A,3),size(F_A,4));  % initialize variable
    %     itc_B_surrogate  = zeros(num_iter,size(F_B,2),size(F_B,3),size(F_B,4));
    %     for i=1:num_iter
    %         randpointA = randsample(Ltimes,Ntrials,1);  % randomly select an onset time point on each trial
    %         randpointB = randsample(Ltimes,Ntrials,1);
    %         for j=1:Ntrials
    %             F_A(j,:,:,:) = F_A(j,:,:,[randpointA(j):end 1:randpointA(j)-1]);  % shuffle
    %             F_B(j,:,:,:) = F_B(j,:,:,[randpointB(j):end 1:randpointB(j)-1]);
    %         end
    %         itc_A_surrogate(i,:,:,:) = squeeze(abs(mean(F_A./abs(F_A),1)));  % compute ITC
    %         itc_B_surrogate(i,:,:,:) = squeeze(abs(mean(F_B./abs(F_B),1)));
    %     end
    %     itc_A.itc = (itc_A.itc - squeeze(mean(itc_A_surrogate,1)))./...
    %         squeeze(std(itc_A_surrogate,1,1));  % Z-score
    %     itc_B.itc = (itc_B.itc - squeeze(mean(itc_B_surrogate,1)))./...
    %         squeeze(std(itc_B_surrogate,1,1));


    % plot
    cfg = [];
    cfg.parameter    = 'itc';
    cfg.showlabels   = 'yes';
    cfg.layout       = 'C:\Users\custs\workspace\matlab\congcong\layout\biosemi64.lay';
    % cfg.zlim         = [-8 8];
    %     figure
    %     ft_multiplotTFR(cfg, itc_A);


    save(['C:\Users\custs\workspace\matlab\congcong\data\mat\' id{subi} '.mat'], 'itc_A','itc_B')
end


%% prepare data for grand-avg and statistical test
close all
clear
cd('C:\Users\custs\workspace\matlab\congcong\data\mat')
id={'103','106','204'}
Ns = length(id);

% prepare the data structure (cell structure)
allsub_A = cell(1,Ns);
allsub_B = cell(1,Ns);

for subi=1:Ns
    load([ id{subi} '.mat']);

    % 将所有被试的数据组织成 cell array
    allsub_A{1,subi} = itc_A;
    allsub_B{1,subi} = itc_B;
end


%% calculate grand average for each condition
cfg = [];
cfg.parameter = 'itc';
cfg.channel   = 'all';
cfg.toilim    = 'all';
cfg.foilim    = 'all';
itc_A_grand = ft_freqgrandaverage(cfg, allsub_A{:});
itc_B_grand = ft_freqgrandaverage(cfg, allsub_B{:});

%%% plot grand average of TFR
cfg = [];
cfg.parameter    = 'itc';
cfg.showlabels   = 'yes';
cfg.layout       = 'C:\Users\custs\workspace\matlab\congcong\layout\biosemi64.lay';
% cfg.zlim         = [-8 8];
figure
ft_multiplotTFR(cfg, itc_A_grand);
figure
ft_multiplotTFR(cfg, itc_B_grand);
% 在界面上选择多个电极，点击后，画出的TFR是这几个电极的平均

%%% plot topography of TFR of effect of A vs. B
% difference TFR(两条件相减)
A_vs_B = itc_A_grand;
A_vs_B.itc = itc_A_grand.itc - itc_B_grand.itc;

cfg = [];
cfg.parameter    = 'itc';
cfg.showlabels   = 'yes';
cfg.layout       = 'C:\Users\custs\workspace\matlab\congcong\layout\biosemi64.lay';
% cfg.zlim         = [-8 8];
figure
ft_multiplotTFR(cfg, itc_A_grand);

% 画出差异TFR的地形图
cfg = [];
cfg.parameter = 'itc';
cfg.comment   = 'no';
cfg.layout    = 'C:\Users\custs\workspace\matlab\congcong\layout\biosemi64.lay';
cfg.xlim      = [0 1.5];   % time limit, in second
cfg.ylim      = [1 30];       % freq limit, in Hz
% cfg.zlim   = [-5 5];
figure; ft_topoplotER(cfg, A_vs_B);




chan2use ={'P3','P4','P5','P6'}; %%'PO4',,'P3','P4','P5','P6',,'PO7','PO8','P7','P8','FZ','FCZ','CZ'
%%%p4 p6 p8 po4 po6 po8 以及对侧效果
% chan2use = {'P7','P8'};
for h=1:length(chan2use)       % find the index of channels to use

    figure %%%这一段有问题，164-168
    cfg.channel      = chan2use(h);
    subplot(2,2,1),ft_singleplotTFR(cfg,  itc_A_grand); title('amci') %%显示单个
    subplot(2,2,2),ft_singleplotTFR(cfg,  itc_B_grand); title('hc')
    subplot(2,2,3),ft_singleplotTFR(cfg,  A_vs_B); title('amci-hc')
    % dlmwrite('PO8_300-450_young.txt',power,'\t')
    time  = itc_A_grand.time;
    freq  =itc_A_grand.freq;
    chan  = itc_A_grand.label;

    % define time window
    %  timewin      = [-0.4 1.2];
    timewin     = [0.15 1.5];  %%[0.3 0.45]      [0.15  0.25]

    timewin_idx  = dsearchn(time', timewin');
    % define frequency window
    % freqwin      = [1 30];  % theta band
    freqwin = [8 13];
    %%%Alpha 波：8–13 Hz；Beta 波：13–30 Hz

    freqwin_idx  = dsearchn(freq', freqwin');
    % define ROI (channels)

    chan2use = {'PO3'};    %%%p4 p6 p8 po4 po6 po8 以及对侧效果
    chan_idx = zeros(1,length(chan2use));
    % find the index of channels to use
    ch = strcmpi(chan2use(h), chan);
    chan_idx = find(ch);

    clear P F
    % extract TFR over these time window, freq window, and ROI, for each condition and each subject
    file_name={'103','106','204'};
    Ns = length(file_name);
    pow = zeros(Ns,2); % initialize variable， 此例中，有Ns个被试，2个条件
    for i=timewin_idx(1):timewin_idx(2)
        for j=freqwin_idx(1):freqwin_idx(2)
            for k=1:Ns
                pow(k,1) = allsub_A{1,k}.itc( chan_idx, j, i );
                pow(k,2) = allsub_B{1,k}.itc( chan_idx, j, i );
                power_h(k,j,i)=allsub_A{1,k}.itc( chan_idx, j, i );
                power_f(k,j,i)=allsub_B{1,k}.itc( chan_idx, j, i );
            end
            [p_value, table] = anova_rm(pow,'off');  %% select all subjects，all conditions，channel Cz for each time point to perform repeated measures ANOVA %%anova_rm会无法识别，需要下anova文件
            P(i-timewin_idx(1)+1,j-freqwin_idx(1)+1)=p_value(1); %% save the p value from ANOVA
            F(i-timewin_idx(1)+1,j-freqwin_idx(1)+1)=table{2,5};%% save the F value from ANOVA
        end
    end

    [m,n]=size(P);%%%  做FDR校正后
    P_ttest1=reshape(P,[m*n 1]);
    P_ttest2= mafdr(P_ttest1,'BHFDR',true);
    P_ttest3=reshape(P_ttest2,[m n]);

    p_fdr=P_ttest3;
    F(find(P>=p_fdr))=0;
    subplot(224), h1=imagesc(timewin,freqwin,F);
    set(gca,'YDir','normal');colorbar;title('p-value with FDR');
    fn2= char(strcat(chan2use(h),'.jpg'));
    saveas(gcf, fn2);
end





%% 导出数据用于进一步统计分析
% 根据已有假设或前一步得到的总平均图，选择时间窗口、频率窗口、电极点（ROI）
time  = itc_A_grand.time;
freq  = itc_A_grand.freq;
chan  = itc_A_grand.label;

% define time window
timewin      = [-0.4 1.2];
timewin_idx  = dsearchn(time', timewin');
% define frequency window
freqwin      = [1 30];  % theta band
freqwin_idx  = dsearchn(freq', freqwin');
% define ROI (channels)
chan2use = {'FZ','FCZ','CZ'};
chan_idx = zeros(1,length(chan2use));
for i=1:length(chan2use)       % find the index of channels to use
    ch = strcmpi(chan2use(i), chan);
    chan_idx(i) = find(ch);
end

% extract mean TFR over these time window, freq window, and ROI, for each condition and each subject
file_name={'103','106','204'} ; %gono3-hc
Ns = length(file_name);
power = zeros(Ns,2); % initialize variable， 此例中，有Ns个被试，2个条件
for subi=1:Ns
    pow1 = allsub_A{1,subi}.itc( chan_idx, ...
        freqwin_idx(1):freqwin_idx(2), timewin_idx(1):timewin_idx(2) );
    pow2 = allsub_B{1,subi}.itc( chan_idx, ...
        freqwin_idx(1):freqwin_idx(2), timewin_idx(1):timewin_idx(2) );
    power(subi,1) = squeeze(mean(mean(mean( pow1  ))));  % 提取第一个条件的数据
    power(subi,2) = squeeze(mean(mean(mean( pow2  ))));  % 提取第二个条件的数据
end
dlmwrite('itc.txt',power,'\t')  % 保存到txt文件中(用excel打开)，用于进一步分析
% 如有其他时频窗口或ROI，也类似操作
% 当选择了多个时频窗口或ROI时，进行了多次比较，此时需要对p值进行校正（eg, FDR）






%%%%%%%%%%%%T检验
% clear P F
% % === 设置参数 ===
% file_name = {'110','119','122'};  % 被试编号（AMCI或HC）
% Ns = length(file_name);
% pow = zeros(Ns,2);  % 存储每个被试两个条件的 ITC 值
%
% % 遍历时间窗与频率窗
% for i = timewin_idx(1):timewin_idx(2)
%     for j = freqwin_idx(1):freqwin_idx(2)
%         for k = 1:Ns
%             pow(k,1) = allsub_A{1,k}.itc(chan_idx, j, i);  % 条件A
%             pow(k,2) = allsub_B{1,k}.itc(chan_idx, j, i);  % 条件B
%             power_h(k,j,i) = pow(k,1);
%             power_f(k,j,i) = pow(k,2);
%         end
%
%         % 使用成对 t 检验替代 ANOVA
%         [~, p_value, ~, stats] = ttest(pow(:,1), pow(:,2));
%
%         % 保存 p 值 和 t 值（代替原来的 F 值）
%         P(i - timewin_idx(1) + 1, j - freqwin_idx(1) + 1) = p_value;
%         F(i - timewin_idx(1) + 1, j - freqwin_idx(1) + 1) = stats.tstat;
%     end
% end
%
% % === FDR 校正 ===
% [m, n] = size(P);
% P_vec = reshape(P, [m*n 1]);
% P_fdr_vec = mafdr(P_vec, 'BHFDR', true);  % FDR校正
% P_fdr = reshape(P_fdr_vec, [m n]);
%
% % 将不显著的 t 值设为 0
% F(P > P_fdr) = 0;
%
% % === 可视化统计结果 ===
% figure;
% imagesc(time(timewin_idx(1):timewin_idx(2)), ...
%         freq(freqwin_idx(1):freqwin_idx(2)), F);
% set(gca,'YDir','normal');
% colorbar;
% xlabel('Time (s)');
% ylabel('Frequency (Hz)');
% title('ITC difference (t-test with FDR)');







